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P-l' 

An algorithm is presented for approximating arbitrary powers of a black box unitary operation, 

■ W*, where t is a real number, and IA is a black box implementing an unknown unitary. The complexity 
£SJ . of this algorithm is calculated in terms of the number of calls to the black box, the errors in the 

approximation, and a certain 'gap' parameter. For general IA and large t, one should apply U a total 
. of \t\ times followed by our procedure for approximating the fractional power U t ~ L*J . An example is 

CIh' also given where for large integers t this method is more efficient than direct application of t copies 

of IA. Further applications and related algorithms are also discussed. 
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I. INTRODUCTION 



, Suppose an n-qubit unitary IA is implemented by evolving (or simulating the evolution of) a time-independent 
Hamiltonian H for a period of time r = 1, that is, U = e~ lH . Then for any t S R+, one can implement IA 1 by simply 



evolving the Hamiltonian for a period of time t. For example, if t = k, then a square root of IA, e l 2 H : could be 



m 

00 . implemented in this way, and in such model of computation the cost would be half of the cost of implementing IA. 

In our work, we question what can be done if U is realized in some other way, such as a non-trivial sequence of 
time-dependent Hamiltonians, or a quantum circuit. In other words, consider the situation when IA is given in the 
form of a black box. The goal is to implement real valued powers t of this unitary operation by making use of the 
, multiple copies of the black box implementing IA. The complexity of such a procedure is measured in terms of the 
\ total number of calls to the unitary. 

It is possible to find the t th power of an unknown unitary by first performing a sufficiently precise complete quantum 
process tomography of a 2" x 2™ dimensional unitary U, which uses 0(4") calls to the unitary with various input 
states 0, Q and measurements to achieve IA with constant precision. The exponential scaling with n is necessary. 
In particular, a lower bound on the number of calls to the unitary can be derived from constructing an e— net over 
unitaries [2^ . If the space of unitary operations is divided up into balls of radius e then the total number of unitaries 
that can be specified up to this resolution is (|) , where c is a constant. To specify each of these requires 4™ log(|) 
bits. To discover information about IA one might supply states to the unitary tensored with the identity operation of 
the same dimension (a larger dimension does not help) and perform measurements on the output states. The states 
will have dimension 4™, and by Holevo's bound one could hope to obtain no more than 2n bits of information for each 
of these calls to IA (in general we will get less than this). Therefore, it is not possible by tomography to use fewer 
than n (£ log(f )) calls to U. Even if we allow for error in some fraction 8 of the basis states, this still cannot reduce 
the required number of calls to any function polynomial in n 0]. 

In this paper, we report on a more efficient approach to implementing powers. In particular, we present an algorithm 
for approximating with error e (using the trace norm) any constant power of an unknown unitary using only 0(j log i) 
calls to the unitary itself, assuming a certain 'gap' parameter g we explain later satisfies g £ 0(e) (otherwise, the 
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complexity is O(Mogi) calls; see Section Til Bl for more detail). Note that this complexity is independent of the 
number of qubits n. 

There has been some relevant prior work concentrating on the relation between discrete and continuous oracles 
for various problems. Farhi and Gutmann Q introduced the concept of (continuous) Hamiltonian black box oracles 
for quantum computing. Ioannou [H], and Roland and Cerf 0] consider the problem of simulating a Hamiltonian for 
Grover search Q on a discrete computer. Mochon Q extends this to a more general setting where he is concerned 
with finding lower bounds for discrete oracle problems by considering them in the Hamiltonian setting and mapping 
them to the problem of finding geodesies in manifolds. In particular, he considers one-item Grover search and oracle 
interrogation, highlighting the case of computing the XOR function on a hidden bit string. In these cases, he is 
able to exploit symmetries in the oracle problems to solve what is otherwise a very difficult problem. Our focus is 
different. Firstly, we do not allow for an oracle to be applied a fractional amount of time (instead, we are concerned 
with simulation of fractional oracles in the case when an oracle is only allowed to be applied an integer number of 
times). Secondly, we present an algorithm, with our focus on (constructive) upper bounds. Thirdly, we do not make 
symmetry assumptions on our oracles. Finally,we do not restrict the eigenvalues of these unitaries to ±1 as in the 
Grover search and XOR case, or the work in [9( which shows how to simulate a general continuous Boolean oracle 
evolving for a total of time T with O (^^j^p) discrete oracle queries. 

The remainder of the paper is organized as follows. In Section [IT] we describe the basic construction for efficient 
computation of the powers of IA, discuss our assumptions and constraints, give the complexity of the computation, and 
calculate the precision of our approximate solution. The full calculation of these results is included in[Bj An example 
of a special case in which t is a large positive integer and our algorithm is more efficient than direct application of 
the t copies of the unitary is given in Section IIIII Related applications of our algorithm to computing fractional 
Fourier transform and noise filtering arc explored in Section IIVI Finally, we discuss the implications of this work in 
the concluding section. 



II. THE ALGORITHM 

In this section, we describe our basic construction. For any unitary IA on a finite dimensional state space, consider 
its spectral decomposition U = PKP\ where P is a unitary matrix composed of the eigenvectors of IA, and A is a 
diagonal matrix containing the eigenvalues of IA. For such a decomposition, powers of IA may be computed as follows: 

U* = PA*P f . (1) 

A key observation is that in the quantum case one does not need to actually implement the basis change operators P 
and P^ in order to exploit the above feature of the spectral decomposition. 
For now, let us assume that t is real number between and 1. 

Our algorithm is described in three stages. In the first stage, approximations to the eigenvalues of IA are calculated 
using an eigenvalue estimation algorithm. Let us initially assume that the eigenvalues are of the form e 2mXk with 
Afc = -|lr for some integer 1^ € {0, 1, . . . , 2 m — 1} (so that we can initially ignore the effect of precision errors). In the 
second stage, phase shifts arc applied to the eigenstates of IA. The third stage uncomputes the eigenvalue estimation. 
Refer to Figure [T] for details. 

Suppose |^) is a quantum state on which we wish to perform the transformation IA 1 . Note that any state |\&) is a 
superposition of eigenvectors, \ipk), oiU, \^f) = ^2 k at IV'fe)- Prepare m ancilla qubits in the "all-zeros" state, where 
to is a user chosen parameter that characterizes the complexity of and errors in our algorithm: its role is discussed in 
more detail in El These ancilla are input to the Quantum Fourier Transform (QFT) circuit [101 ] (or equivalently, since 
applied to the all-zeros state, a circuit consisting of only Hadamard operations on each qubit) and later they are used 
as the controls for controlled-W operators, while the state \^) is the target. The notation c — W denotes the operation 
\j) \4>) h- ► \j)W Using the eigenstate expansion of yS>) in the target, it can be shown that such a controlled 
operation causes a phase kick-back to the control register, so that the total state becomes ^ fc J^j afce 2 ™^ |j) \ipk)- 
Next, an inverse QFT operation is performed on the ancilla register. The result left on the ancilla register are 
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FIG. 1: The quantum circuit diagram outlining the implementation of the algorithm that raises an unknown unitary U to some 
power t. Stage I performs eigenvalue estimation, where value I in the first register corresponds to the eigenvalue parameter 
estimate A = l/2 m . Stage II then uses this eigenvalue estimate to induce a corresponding phase shift e 27Tdt / 2 . Stage III 
uncomputes the eigenvalue estimation step. In the case that the eigenvalue parameters Xk are of the form Z^/2 m , the eigenvalue 
estimation is exact and the uncomputation step returns the control register back to the initial all-zeroes state. In general, when 
we do not have this assumption on the form of the eigenvalues, we perform a slightly more complicated eigenvalue estimation 
algorithm that gives an estimate with error at most l/2 m with probability in 1 — 0(l/2 m ). 



estimates of the eigenvalue parameter, Afe, specifically l k = 2 m Afc, in superposition (as dictated by the input vector 
| VP)). For a more detailed description of the eigenvalue estimation algorithm see 0. 

For stage two, we apply the controlled operation c — 7£(27r£fci/2 m ) to induce the phase shift e 27 ™^*/ 2 ™ when the 
eigenvalue parameter estimate £k is in the control register. This computation corresponds to the exponentiation of 
the diagonal matrix in the spectral decomposition formula. There are a few natural ways to implement this step, as 
shown in 



11( or (12 1 . 



Finally, the state of the control register (containing the £k values) is uncomputed back to the all-zeros (input ancilla) 
state. This completes the third and final stage of the computation. 

This leaves the registers in the final state |0) CS> ^2 k ctke 2mXkt \ipk): which is simply the result of the application 
of U l to the state | \P) . When we drop the assumption that the eigenvalue parameters Afc are exact ratios of some 
integers i k divided by 2 m , we will in general get precision errors in the eigenvalue estimation. This imprecision in the 
resolution of X k means the phases applied to the state are not exact, which further implies that the uncomputation 
step does not return the ancilla register precisely to all zeroes. Still, by a proper choice of parameters, and a "gap" 
assumption discussed below, the errors can be managed — in particular, we prove that they are exponentially small 
in the parameter m. This error reduction is done by choosing an eigenvalue estimation algorithm that outputs an 
estimate with error at most with probabilityjn 1 — O(i). This can be done with 0(m) repetitions of the standard 
QFT-bascd eigenvalue estimation algorithm 10( and applying the Chernoff bound. This uses 0(m2 m ) calls to the 
black box. Relevant calculations may be found inlBl 



A. Underlying Assumptions 

Our algorithm, as formulated in the previous section, asserts that in addition to the black box implementation of hi 
itself we also have access to a black box implementing controllcd-W and a black box implementing U^ 1 . Furthermore, 
we also must assume, as discussed in more detail in the following section, that spectrum of matrix hi has a small 
gap, in particular, that U has no eigenvalue A = e 1 ^, where <j> e (27r(l — g), 2tt) for some value g. For simplicity, we 
assume that g > otherwise, we would need to replace some of the 0(2 m ) terms with 0(max{2 m , |}) (the reason 
is that we need to do eigenvalue estimation with precision less than the gap size with high probability). In this case, 
we achieve an approximation with error in O(^). Note that there are no other assumptions about gaps elsewhere in 
the spectrum. Furthermore, an assumption of this form is essentially necessary in the case of worst-case complexity 
(as we explain in Section [II B| . which is what we are interested in. Average-case performance is discussed more inlDl 

Not all of the above restrictions are necessary, but they facilitate a clear and transparent description of the algorithm 
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and its analysis. In some cases of interest, the individual restrictions may be dropped. In particular, in cases where only 
the black box implementation of lA is given (and not the controlled-ZY) , one can exploit the technique for implementing 
a controlled-e -1 ^?/, where e 1 ^ is a random eigenvalue of lA, and use it in the algorithm to find arbitrary powers, as 
outlined in 13[ (see for details). In some special cases, we can remove the assumption of having a black box for 
lA^ 1 . For example, for non-integer values of t > 2 m , we can approximate with error in O(j) using 0(t) calls to the 
controllcd-W, without any use of a black box for lA -1 (see[C]). We note that the above restriction on the form of the 
spectrum is naturally satisfied in many practical unitarics, including the quantum Fourier transform, standard oracles 
\j) i— > (— \j) for computing Boolean functions /(•), and oracles such as those implicitly used in the adiabatic 
algorithms of Farhi et al. |14[ that reversibly compute a function with a finite discrete spectrum. 

It is also worth noting that, for non-integer values of the power t, there are many operators that one might reasonably 
call the t th power of lA. In other words, for any given lA, there are many Hamiltonians H such that lA = e~' lH , and 
thus one could naturally define lA 1 as e~ %Ht for any one of these Hamiltonians. For instance, consider t = \ and the 
identity matrix I = ( Q 1 ). Let us find all unitary matrices B such that BB = I. In the most natural understanding, 

each such B is a possible square root of I. Firstly, matrices I and ^ 1 _ < | ) ^) are both square roots of I. This happens 
because there are two possible square roots of the complex number 1, which is the eigenvalue of I. In addition, each 
eigenspace of a unitary matrix may be broken into subspaces such that different square roots of the corresponding 
eigenvalue may be used to construct square roots of a given matrix. In our case, this helps to construct two more 
square roots of the form ± ( n _j. ). Furthermore, every matrix of the form ± ( c ° s b a . e ^ sma \ anc j its transpose, 



— 1 J ' ' J \e b sin a —cos a 

where a and b are real valued parameters is a square root of the operator I. For the choice of parameters a = ^ and 
6=0 this allows to construct root ( ° J) , also known as the Pauli-X matrix. Roots of the identity matrix in higher 
dimensions get more complicated. Thus, the problem of finding fractional powers is not completely straightforward 
to define. In our solution, we construct fractional powers based on the spectral decomposition, choosing the primitive 
complex fractional power of every eigenvalue (though, in principle our method allows us to break the eigenspaecs into 
subspaces and take different roots of the eigenvalue for each subspace). Thus, in a sense, the fractional power we 
construct is the most natural one, and also, it is the fractional power of a matrix most commonly defined in linear 
algebra textbooks. With the above example, our algorithm implements the I root of I, though modifications are 
possible that would allow us to explore a wider range of square roots of I. 



B. Optimality 

For such a general black box and assuming a spectral gap of size g = ^r, we cannot compute IA^ with constant error 
with fewer than £l(2 m ) calls to lA. To see why this is the case, suppose we are given either lA = I or lA — lA\-\/2^ = 
|0) (0| + e~ 27 ™/ 2 |1) (1|. It follows from the following simple lemma that Q(2 m ) calls to lA are necessary to correctly 
guess which lA we were given with probability greater than |. 

Lemma 1: Suppose we are given a black box lA that implements either I or lA\~s = |0) (0| + e |1) for some 
small 5 > 0. Any algorithm that correctly guesses the identity ofU with probability at least | ; {or any prior distribution 
of the two possible values of 14, must make evaluation s of lA. 



This can be easily proved, e.g., by the "hybrid" method la ]. 

Note that the root ^_ 1/2m (|0) + |1}) = |0) - e^ 7 ™/ 2 '" |1) is almost perfectly distinguishable from li (|0) + |1)) = 
|0) + Thus, the power finding algorithm must take 0(2 m ) calls to lA in order to compute IA^ with (small) constant 
error. The same holds for any power t € (0, 1) that is bounded away from from both and 1. 

Thus, while more efficient ways to compute fractional powers of lA are possible in special cases, in general one 
cannot do much better than what we have described. In the case when < .g < 7r, we can still obtain a lower bound 
of f2(2 m ) queries by letting lA be either Z or |0) (0| — e~ 27 ™/ 2 |1) (1|, and computing square roots of lA 2 as a means 
of distinguishing the two cases. 

This same argument also shows why a gap assumption is necessary. Without the gap assumption, there is no upper 
bound on the worst-case complexity of computing the square root of a black box unitary, since such a black box could 
be used to distinguish I from IA1-5 = |0) (0| + e~ i<5 11) (1| for arbitarily small values of S. As pointed out earlier, if 
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we have a known gap g without a promise that g £ f2(e) = f2(l/2 m ), then a complexity of O(Mog-j-) suffices, and 
the example above shows that fi(-) queries are necessary. If g is unknown then there are several cases one might 
consider. If we are promised a lower bound < <7 m in < 9-, then in the worst case, we must use 0{-^— log-i) calls 
to the unitary. If instead we have an means of reliably recognizing when the algorithm has succeeded, then trying 
a sequence of lower bounds that decreases exponentially until the algorithm succeeds will also give a complexity in 
0{ min | g e } log \), which is the same complexity as when g is known. If, however, we do not know g and have no means 
of recognizing when the algorithm has succeeded, then no upper bound is possible. 



III. LARGE POWERS AND AN EXAMPLE OF EXPONENTIAL IMPROVEMENT 



When it comes to computing large powers of hi by this method, precision in the eigenvalue estimates becomes a 
limiting factor. If instead of the exact eigenvalue e 2m Afc we have its estimate e 27 ™ Afc , where | Afe — Afe | = e, then mapping 
\ipk) l— ► e 2mXkt will lead to an error in the phase parameter of size te. For arbitrary hi, it takes 0(-) calls to 
hi in order to get an estimate with error at most e with high probability. Thus to achieve an error of 5 — te in the 
application of W ! , we need e = 6/t precision in the estimate, amounting to 0(|) calls to hi. For t > 1, a better thing 
to do in this general case is to first apply hi a total of \t\ times, and then use the method from the previous section 
in order to compute the remaining fractional power of hi. 

However, if for special cases of hi we are able to get smaller errors in the eigenvalue estimates, we can do much better. 
For example, for black boxes that map \j) i— ► ( — l) Xj \j) for some j £ {0, 1}™ and X1X2 . . . Xj . . . Xn £ {0, 1} , it is 
easy to exactly determine the eigenvalue parameter Xj with only one call to the black box. As discussed later, it is 
easy to exactly determine the eigenvalue of an eigenvector of the quantum Fourier transform, since the QFT has order 
four. In both of these examples, our technique does not help with computing higher powers of hi since hi has very 
low order, say r, and thus computing hi 1 is equivalent to computing hi 1 mod r , where t mod r is the unique element 
of t — rZ that lies in in the interval [0,r). We next construct an example where hi has very high order, but one can 
still compute precise estimates of its eigenvalues with relatively few calls to hi, and thus exponentiate hi with much 
fewer than t calls, even for t less than the order of hi. 

Suppose the eigenvalues of the desired unitary 14 are all of the form e p *= , for k £ {1, 2, . . . , b\, where tk is an integer, 
Pk £ {pi = 2,£>2 = 3,^3 = 5, ...,pb\pi = i th prime}, and b is a natural number. Then hi B = I, where B = 
and (assuming that each pj» occurs at least once with a non-zero i&) no lower power of hi is equal to the identity. If t 
is on the order of B, then the number of calls to the unitary required for the straightforward implementation of U 1 is 
in 0(B) = O (U b k=1 k log k). This product is called a primorial 0, and it is in 0( e ( 1 +°( 1 )) f ' 1 °g b ) = 0(M 1+o(1 » b ). 

We next demonstrate how our basic algorithm introduced in Section |TT] may be updated to approximate hi 1 with 
exponentially fewer calls to hi than the number of calls it takes to obtain hi 1 via a direct application of t unitaries. 

The method used for this problem employs the continued fraction algorithm, which was also used in the order- 
finding part of Shor's algorithm [17]. The eigenvalue estimation algorithm will return an estimate t^ft of the eigenvalue 
parameter to within an error at most -^r with probability at least ^ . The error probability can be amplified down 
to by repeating the estimation 0(m) times and taking the majority vote. If m is chosen such that m £ 0(\ogb), so 
that 2 m > 2pbPb-i, then the fraction ^ will be the only fraction, or convergent, in the continued fraction expansion 
of that has distance at most from ^ and has a denominator at most pb- 



h-l h h_ h+1 

FIG. 2: A diagram illustrating the procedure. The estimate is ^r. The smallest gap between any two eigenvalues is j-and 
we require this to be larger than ^ir to ensure that only one rational fraction with denominator no greater than pt, will be 
found in the region to . 
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This convergent ^ can be found with 0(m) arithmetic operations over integer numbers with at most logm bits 
each (e.g., no more than 0(m 3 ) binary operations using trivial algorithms for division and multiplication) using the 
continued fraction algorithm. Once we have the exact value of the eigenvalue parameter with probability ^r, 
we can correctly compute any positive power t of U with error in O(^) using poly (b, log t) calls to U and other 

elementary operations. This is done by mapping, in superposition, \£k) \pk) \ipk) l— * e 27ri * p fc + °' 2 ™ •* \£k) \pk) \ipk), which 
can be done by standard reversible computing methods and m-bit arithmetic operations (the O(^) error is due to 
precision in the arithmetic operations). Then, as usual, one uncomputes the computation of \pk). As alluded to 
before, this gives an overall complexity of 0(m2 m ) = 0(b\ogb) to perform the operation up to error O(t^). This 
follows from (|B11|) in [B] 

This complexity is polynomial in logi (vs poly(t)) because we only need to compute an m-bit approximation to 
t£k/pk- This also allows us to uncompute the eigenvalue estimates found in the first stage of the algorithm without 
ever needing to use the operation . Note that t could be very large, for example t ~ B/2 = 2 b ~ 1 . In this case, we 
can apply U* to any state of the appropriate dimension with the output state containing an error term with probability 
amplitude 0(\) with only 0(6 log 6) <C t G 0{2 b ) calls. The error probability can efficiently be made exponentially 
small by repeating the phase estimation part. 



IV. FRACTIONAL QUANTUM FOURIER TRANSFORM 

In classical computation, roots of the Fourier transform are used for efficient filtering of the frequency noise that is 



some non-constant function of a conjugate variable, such as time [181 ] . Given the role the QFT plays in the quantum 
computation, and the importance of roots of Fourier transforms in classical computation, there may as well be other 
interesting techniques based on the fractional powers of the quantum Fourier transform. 

Note that QFT" 1 = I. Let us illustrate the computation of a fractional power t, t < 1 of the QFT using the algorithm 
introduced in Section[n] First, note that the QFT has four distinct eigenvalues: +1, i, —1 and —i. This means that 
we need only two ancilla bits for the eigenvalue estimation part of the exponentiation algorithm, and such estimation 
results in the exact values. During the eigenvalue estimation algorithm, the controlled-QFT is applied three times 
(2° = 1 with the first ancilla qubit as a control, and 2 1 = 2 more with the second ancilla as a control). We need three 
more applications of the controlled-QFT in the third stage of the algorithm to uncompute ancilla. This means that 
the query complexity of the computation of a fractional power of the QFT is exactly six queries to the controlled-QFT. 
Moreover, during such computation we introduced no errors due to the approximations of eigenvalues, and thus, given 
all gates are perfect, the computed fractional power is exact. 

This specific example of fractional QFT was alternatively derived as a special case of the work of Klappencckcr 
and Rotteler [l^ ] , where they address the different problem of how to simulate a known unitary matrix A given other 
operators that generate an algebra containing A. 



V. CONCLUSIONS AND DISCUSSION 

In this work we have further developed the connection between continuous time black boxes (i.e., "Hamiltonian 
oracles") and their discrete counterparts, unitary black boxes. Previous work has mostly focused on important special 
cases of this relationship. Our work pushes the boundary in a new direction, where the unknown hi is a unitary acting 
on an arbitrarily large Hilbcrt space. It also demonstrates constructively that an algorithm with a continuous oracle 
can also be implemented with a discrete oracle that implements the continuous oracle for some fixed time (with some 
assumptions on the Hamiltonian encoded by the oracle; e.g., if the energies, the eigenvalues of the Hamiltonian, are 
between and E max , then the fixed time interval should be bounded below 2ir/E max so that our gap assumption is 
satisfied). 

For the case that t is a positive non-integer, the complexity of this algorithm is substantially higher than that in 
the continuous setting, where we have access to a Hamiltonian H satisfying hi = e~ lH . If we have direct access to 
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the Hamiltonian, H , then in principle we can easily perform any unitary of the form e~ lHt with a cost t in terms of 
calls to the Hamiltonian. For an algorithm with total query complexity t < 1 in the Hamiltonian black box model, it 
is hard to imagine a reasonable way to effect e~ lHt with fewer than one evaluation of U, and indeed, our algorithm 
makes O(-Mogi) calls to U to effect e~ lHt with precision e (wc show a corresponding lower bound in section Iff Bp . 
which is notably lower than the number of calls required in the obvious approach of doing process tomography on U, 
approximating U l , and then synthesizing a circuit to implement such an approximation. In particular, our approach 
docs not depend on the dimension of the Hilbcrt space that U acts on. 

We have presented an algorithm for finding the powers of unknown unitary operations raised to any real power t > 0. 
We found the complexity of the algorithm to be in 0( [t\ + 2 m m) where the error is in 0( ^r). For large integers t we 
presented a non-trivial example where our method is exponentially more efficient than the direct repeated application 
of U a total of t times. Note that the same exponential speed-up is possible also in the continuous case if the 
Hamiltonian has some structure that enables eigenvalue estimation much better than the standard worst-case limits; 
in such cases, our procedure, using a modest amount of ancilla workspace, can simulate the evolution of H for a time 
r in time polynomial in log(r). 

Our algorithm could be of practical use in a situation where we have very precise clocks and control, say with 
error less than a very small positive e, but for some technical reason wc must run the unknown Hamiltonian for a 
minimum time T m j„ >> e (e.g., because of a minimum reset time on some critical piece of apparatus, like a detector, 
used to simulate or execute the Hamiltonian), and T m i n is still smaller than 2-K/E max (assuming the energies are 
non-negative). In this case, there is a smallest amount of time for which a given Hamiltonian can naturally be run. 
If we wish to run the Hamiltonian for less than that amount of time, then our method could be used. 

Further, if the Hamiltonian that solves some problem is time dependent, as is often the case in quantum circuits or 
in models based on the permutation or braiding of particles on a surface, it may be of interest to apply some fraction 
of this total transformation to turn the solution to the problem into a subroutine of a larger algorithm. The method 
we propose provides a general way of obtaining this fractional application. 

Our technique may also be applied to compute other functions of U, by replacing each eigenvalue e lXk with e 4 -^ Afc - ) 



for any function /(■). This may be useful in contexts such as Childs' method [20( for simulating continuous time 
quantum walk Hamiltonians with discrete time walks. This requires applying phases that are the sines or arcsines 
of the eigenvalues of the original Hamiltonian. Our method could achieve this for an unknown Hamiltonian to be 
simulated. The properties of the function /(■) will determine the efficiency as a function of the precision. In the case 
when /(•) is continuous with bounded first derivative, and 27r-periodic, then we can also drop the assumption that 
there must be a gap in the spectrum. 

In addition, there arc a number of classical cryptography schemes that rely on the computational difficulty of finding 
roots modulo some number, such as RSA. One might naturally attempt to construct quantum cryptography schemes 
where knowledge of the "secret key" corresponds to the ability to compute a fractional power (i.e., a root) of some 
unitary operator. This algorithm provides a way of implementing such roots, implying serious limitations for such an 
approach in quantum cryptographic schemes. 
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APPENDIX A: CONTROLLED UNITARIES 



Note that the eigenvalue estimation algorithm requires the use of a controlled unitary c — hi. If we are given 
an eigenvector of hi with eigenvalue +1, then this transformation can be implemented using a controllcd-SWAP 
(equivalently, a series of Fredkin gates 21|), as was shown by Kitaev 1 1 31 ] and is illustrated in the following circuit 
diagram: 



1 91/ control 



92 



SWAP 



It is easy to verify that the above circuit implements 



|0) 



control I ^2 /target 



1^ 



u 



SWAP 



|0) 



control I ^/ target 



lb 



control 



I 92, 



target 



lb 



control 



1 92) target) 1^* 



In other words, up to a global phase, this implements the controlled- (e~ l7VlXk U). If \k = 0, i.e., \tpk) has eigenvalue 
+1, then we have implemented the c — hi. Otherwise, we have a very similar operation, which may or may not be 
sufficient, depending on the application. If, for example, any c — (e^hi) will suffice, as long as we use the same one 
consistently, then this is easy to achieve. 

We do not need to assume that we know the eigenvectors {iV'fc)} of hi or are given a specific eigenvector. Instead, 
we can use any state in the target register. For simplicity, we can consider the completely mixed state, which can be 
decomposed as an equal mixture over the eigenstates. Provided we keep the state that is in the target register, and re- 
use the resulting state every time we construct a controlled-^ between these two registers, the phase (though random) 
will be the same for the entire computation. That is, this is equivalent to implementing a c — (e l ^hi) consistently 
throughout the computation, where <p is an unknown random value (whose distribution depends on the weight of the 
cigenspace of e~ l< ^ in the initial mixed state). In essence, the target system serves as a phase reference. 



APPENDIX B: COMPLEXITY AND ERROR ANALYSIS 

We will describe upper bounds on the precision of our algorithm, which is constrained by the lack of perfect precision 
in computing the eigenvalue estimates. We give the analysis for t = i, however the same final bound applies for any 
t g (0,1). 

Let us first describe an "ideal" algorithm that exactly computes the square root of hi. For the ideal algorithm, some 
of the steps will be unnecessary, but we include them to make the comparison with the actual algorithm simple. 
Stage 1 computes a near-optimal eigenvalue estimation (in the sense of nearly optimizing the chance of obtaining an 



estimate with error at most ^r)- We do this by repeating the "standard" eigenvalue estimation algorithm [ll| a total 
of r E 0(m) times and take the majority answer, with the constant chosen so that each eigenvalue estimate e 2nlXk 
of e 2mXk satisfies |e 27 " Afc — e 27rlAfc | < with probability in 1 — 0(^^). One can fine-tune the optimal eigenvalue 
estimation procedure further (e.g., [22|), however, this simple procedure is sufficient. 
This eigenvalue estimation procedure is described by the transformation 

io> y, a * 1^) ~ E a A E iy> i A y) J ( B1 ) 

k k \ y / 



9 



where the values y £ 7L r 2m (where TL^^ = {0, 1, . . . , 2 m — 1}) arc r-tuplcs of integers, and A y is the value obtained by 
taking the most commonly occurring integer j/ mo d e i n the r-tuple of integers y (with some convention for breaking 
ties) and letting A y = y m odc/2"\ 

In the ideal Stage II, we exactly implement the map 

" E a * (E Py* \y) l A y>) e ^ lXk ' 2 1^) ■ ( B2 ) 

Then Stage III uncomputcs the eigenvalue estimation exactly (since with the ideal Stage II, there is no coupling of 
the control registers with the target registers), leaving us with the ideal output 

^« fe |00...0)e 2 " Afc / 2 |^). (B3) 

fc 

The only difference between this ideal algorithm and our actual algorithm is that Stage II actually effects the 
unitary 

E a « ( E iy> wlw-E^E iy> I A y> e2mXy/2 w*) ■ m 

k V y / k y 

Since the probability that |e 2, ™ Ay / 2 — e 27rlA fc/ 2 | > _L i s m O(^r) (note that here we use the "gap" assumption), we 
can write the state in Stage II as 

i $ > = E^E^iy)i A y) e2 " Ay/2 i^) 

k yes 

+E«*E \y) i A y) e2 ™ Ay/2 1^) • ( B5 ) 

k y£S' 

where the values of y £ S include the "good" values of y that produce estimates A y satisfying |e 27riAy — e 27riA k | < J_ j 
and the "bad" values, y £ S', are the rest. The norm of the bad part of the state is in 0(^^) and the norm of the 
"good" part is 1 - O(^L). 

Let A denote the ideal phase shift operator, and A denote the actual operator. 

We thus have 

ai$)-ai$) = E^E^iy)i A y)( e2 " Afc/2 - e2 " Ay/2 )i^) 

k yes 

E Py* ly> IV) ( e2 " Afc/2 - e2 " Ay/2 ) i^> ■ 

fc yes' 

Thus we can bound the norm of this difference by 

||A |$) - A |$) || 2 < J2 E \ a kPyM 2 \Sk, y \ 2 + E E \<x k f3y,k\ 2 \h, y \ 2 , (B6) 

yGS fc yeS' 

where |<5 fe , y | 2 = |e iAfe / 2 - e iA "/ 2 | 2 . 

We know that EyesEfc K/?y,fc| 2 G 1 - O(^), £ y£S , Efc K/3y,fc| 2 G O(^), and that |4, y | 2 G O(^r) for the 
good values, and at most 1 for the bad values. 

Thus 

||A|d>) - A|<&) || 2 e O (^r) - (B7) 
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Noticing that 

II |«> - \v) ||i = ((u\ (v\)(\u) \v)) > 2(1 - \{u\v)\) (B8) 

and using the equality 

|| \u) (u\ - \v) (v\ || Tr = 2y/l -|(4>| 2 , (B9) 

which implies that 

lll«)(«|-|«)Hlk<2|||«>-|«)|| ai (BIO) 

so ||A(|$) ($|)A — A(|$) ($|)A||tc € 0(2^-). The trace norm of the difference between the two unitary superoperators 
is simply given by taking the maximization over states |$). Since this holds for any state |$), we have 

HA-AlkeO^-L). (Bll) 

This is the case after Stage II of the algorithm. Note that Stage III is a unitary operation that is the same for both 
the ideal and the actual case. The trace norm is invariant under unitary transformations, so the trace norm of the 
unitary operators for the full ideal and actual algorithms have the same bound. In the case of unitary operators this 
is equivalent to the diamond norm. 



APPENDIX C: INVERSES 

Our algorithm for generating U l uses U~ x when the eigenvalue estimates need to be uncomputed. More specifically, 
we start with |0) \tp k ), and we apply (QFT~2 <E> I)c - U j (QFT 2m ® I) to compute 

\ QFT 2m I j) ) IV'fe) ■ (CI) 



j 



We then use the value in the first register to control the phase shift corresponding to \ip k ). For simplicity, let us 
assume that we do not bother to induce a phase shift, but we are still left with the task of uncomputing the "damage" 
we have done. It is trivial if we have access to c — U~ 1 (which can be implemented with ac-W and U^ 1 ), since we 
just apply (QFT 2 J ® I)c - U- J {QFT 2 ™ ® I) to get back |0) \ip k ). 

However, in some cases, we do not need the U~ x operations. We could attempt to circumvent the need for U~ x 
operations by applying (QFT 2m ® I)c - W(QFT 2 J <g> I) to yield 

{QFT 2m <g> I)c - W(QFT 2 J ® 1){QFT 2 2 ® l)c - W (QFT 2 ™ <g> I) |0) \i/> k ) 

= (QFT 2 m (g) I)c - W{M ® l)c - U j (QFT 2 m ® I) |0) \tp k ) (C2) 

where M = QFT 2 J, which implies M\x} = \2 rn - x mod 2 rn ). 

We can easily verify that (M ® I)c — W(M ® I) = c — U 2 " 1 ^^ mod 2 "\ and thus the above state equals 

= (QFT 2m <g> I)(M g> I)c - W 2m " j mod 2 ™(c - U j (QFT 2 m ® I) |0) |^ fc ) 

= (QFT 2 J ® I)c - W 2 " 1 (QFT 2m ® I) |0) |^ fc ) . (C3) 

Note that since the c — U 2 no longer depends on the value of the first register, we can commute it through the 
QFT's to get 

c - U 2m (QFT 2 J <g> I)(QFT 2m ® I) |0) \tp k ) 

= c - U 2 " 1 |0) \^ k ) = e 2 ™ 2 '" A * |0) |V fc ) . (C4) 
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This phase factor of e 27 ™ 2 Afc is in general a problem when the target register is in a superposition of cigcnstates, 
since different eigenstates will pick up a different phase factor. 

In some cases, this is not a problem. One observation, is that such a transformation is equivalent to applying U 2 . 
Thus, if our goal is to apply U l for t = [t\ + 5, < 5 < 1, where t > 2 m , we can apply our algorithm for implementing 
U & , but instead with the above modification. This yields an approximation to U 2m+S with error in O(i). We can 
then apply W*~ 2 "* to complete the approximation of U l with error in 0{\) C 0{-^). 

Also note that, if A& = ^ for an integer i^, then the phase factors e 2T ™ 2 '" Afc = 1 all equal 1, and thus pose no 
problem. 

This problem can also be remedied in any other cases where Afc can be determined exactly (or more significantly 
precisely than error jw), since before uncomputing A&, one can add an additional phase shift of e -27 ™ 2 Xk conditional 
on the value of A&. This will eliminate the final unwanted phase factor of e 2Tri2 Afc associated with the eigenvector 
|^fc>. 



APPENDIX D: COPING WITH NO GAP ASSUMPTION 

If one is only interested in the average-case performance (e.g., averaging over all possible input states to 
according to the Haar measure) , one is still faced with the potential problem that most or all of the spectrum of U is 
in the difficult region (see section Hi B[) of f» 1, 7T < <f> < 2n. 

If one is content with obtaining a square root of e l8 U for a random 8 £ [0,2tt), then one can show that the 
average-case error (for any input distribution) is still in O(^r) for an algorithm that only uses the square root of e %B U 
once. 

For an algorithm that uses (e te U)^ a total of k times, the error could get mag nified to O(Jfc). Thus, one could 
select m' = rn + \2 log k] for each simulation of (e ie U)^ . 

In order to see why we have a k 2 factor instead of k, we present an algorithm where substituting an exact square 
root with our approximate square root yields an error in Q(M^). The value 9 is of the form 2n£/2 m — e, for some integer 
I and non-negativevalue e < l/2 m . Consider an operator U that has every2 m th root of unity e 2 ™^ 2 as an eigenvalue 
with equal multiplicities. Consider analgorithm that does eigenvalue estimation on the identity operator withprecision 
in 0(l/2 m ), and uses quantum searching to find aneigenvector with eigenvalue — 1. Clearly, this algorithm should 
not find such an eigenvalue since all the eigenvalues are +1. Now let U\ — (e ie U)~^ and let U2 (for simplicity, let 
us assume it is unitary, even though in practice it will be a map that is almost unitary) be our approximation to 
(e l9 U)^ for the same 9. Note that Ui(e %e U)^ = I, but U1U2 will have most eigenvalues near +1, and proportion 
1/2™ of the eigenvalues near — 1. Thus, a generalized version of quantum searching (sec [Ejl willgcnerate these (close 
to) —1 eigenvectors with an amplitude in0(^=) using only k evaluations of U\U2- Recall that with the ideal (e l6 lA)^ 
(insteadof U2), there would be no such eigenvectors for the search algorithm tofind, thus the part of the wave function 
that found these —1 eigenvalucsis an error. This implies an error with probability in8(^-). 



APPENDIX E: A GENERALIZATION OF QUANTUM SEARCH 

Suppose we are given a black box that flags an unknown quantum state \<f>) £ Hjy. In other words, \4>) \b) = 
>} \b © 1) and \(f>'} \b) = \<f>') \b) for any \<f>') orthogonal to For brevity, let us directly use the phase flip version 



of this black box \4>) = — \<f>) and \<jjf) = \<j>') for (4>\ 0') = (see Section 8.1 of [10( for a more detailed discussion 
of the relationship between these black boxes) . 

Note that if we have a unitary operation A that maps |00...0) 1— > sin(0) \<j>) + cos(9)\cf>') where (4>\<fi') = 0, 
then we can define the slight generalization of a quantum search iterate, Q = — AC/00. ..oA~ 1 O l f > , where U00...0 = 
I - 2 1 00 . . . 0) (00 . . . 0| and note that Q k |00 . . . 0) = sin((2fc + 1)0) \<f>) + cos((2fc + 1)0) \<f>') . Thus, if we have an 
algorithm A that "guesses" \4>) with probability p = sin 2 (0), then after 0(-^=) calls to Q, we obtain |0) with probability 
very close to 1. 
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In traditional quantum searching, the state \(j>) is known to be a computational basis state, and thus it is very easy 
to produce a unitary A such that | (0| A |00 . . . 0) | = 7=, in particular, any operator that maps 1 00 . . . 0) i— > \ x ) 
will do. 

However, for an arbitrary quantum state |0), it is non-trivial to produce such an initial state. One cannot simply 
use a maximally mixed state, since amplitude amplification requires the inverse of the operation A that generates 
the "guess". We could resort to techniques for generating pure states that for practical purposes are "random" (e.g., 
[23)]), however, for this particular application, it is sufficient to use the following idea. 

Let A act on an N 2 dimensional Hubert space and maps 1 00 . . . 0) 1 00 . . . 0) 1— > J2 X ^ 1*^' Note that for any 
basis {|0,)} 7 this maximally mixed state can be rewritten as Ylj Ty^f \4>j) 14^) (where the coefficients of |0*) are the 
complex conjugates of those of \4>)). Notice that the maximally entangled state has equal support for every eigenvector 
in the basis of eigenvalues of the black box unitary, which we require for the procedure outlined in this section to 
succeed. Let us assume, without loss of generality, that |0o) = \4>)- If we a Pply amplitude amplification, using 
Q = —AUoo...o.o...oA~ 1 (0 ( f > Cg> /), then it is easy to verify that 



Q fe |00...0)|00...0) = sin((2fc + l)0)|0 o ) |0£) 



cos^fc + l^l^ 1 -^^,)^*)], 



(El) 



where sin(#) = -7=. Thus, by choosing k w ^ € 0(--j=), we obtain |0 O ) |0o) w i tn n ^S n fidelity. For applications where 
the extra |0*) system will pose a problem, this approach will not work, and techniques for generating pure states that 
behave like random quantum states could be used instead. 

This technique also works if the black box flags a subspacc of dimension d. Again, assume without loss of generality 
that this subspace is spanned by |0o) ; |0i) ; ■ ■ ■ ; \4>d— l)- Then the maximally entangled state can be rewritten as 




sin(0) V -= |0,) |0*) + cos(0) V -=== |0,) |0*) (E2) 



where sin 2 (6') = jj. 

Thus Q k |00 . . . 0) |00 . . .0) « J2 d jZo 7s l^i) |^) forA,- f yjf , assuming <J jj is much less than one. 



Furthermore, if we applied amplitude estimation 2J| in this case, we would be approximating the value of 9 
arcsin(w -^-), which would allow us to approximate the dimension of the subspace, d. 
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